DESY 96-057 
April 1996 



ISSN 0418-9833 



LEPTOft 6.5 — A Monte Carlo Generator 
for Deep Inelastic Lepton-Nucleon Scattering 

G. Ingelman a ' 6 , A. Edin fo , J. Rathsman 6 

ingelman@desy.de rathsman@tsl.uu.se edin@tsl.uu.se 

a Deutsches Elektronen-Synchrotron DESY, Notkestrasse 85, D-22603 Hamburg, FRG 
b Dept. of Radiation Sciences, Uppsala University, Box 535, S-751 21 Uppsala, Sweden 

Abstract. Physics and programming aspects are discussed for a Fortran 
77 Monte Carlo program to simulate complete events in deep inelastic lepton- 
nucleon scattering. The parton level interaction is based on the standard 
model electroweak cross sections, which are fully implemented in leading or- 
der for any lepton of arbitrary polarization, and different parametrizations 
of parton density functions can be used. First order QCD matrix elements 
for gluon radiation and boson-gluon fusion are implemented and higher order 
QCD radiation is treated using parton showers. Hadronization is performed 
using the Lund string model, implemented in Jetset/Pythia. Rapidity gap 
events are generated through a model based on soft colour interactions. 

1 Introduction 

Deep inelastic lepton-nucleon scattering |l] has played an important role for probing the 
structure of the proton and understanding the fundamental electromagnetic, weak and 
strong interactions at the level of quarks and leptons. With the order of magnitude 
increase in the centre-of-mass energy now available through ep collisions in HERA, this 
line of research will continue to be at the forefront. 

The basic lepton-quark scattering processes have well-defined cross section formulae 
within the electroweak standard model p[. With the inclusion of parton density distri- 
butions and perturbative QCD corrections the problem of practical evaluations become 
quite complex and analytical calculations are only possible in simplified cases or through 
approximations. Normal numerical methods are in many cases possible, but Monte Carlo 
simulation is often preferable because of its generality and applicability to complex prob- 
lems. In the case of the multiparticle hadronic final state, the only viable alternative is, 
in fact, the Monte Carlo method. 

The present program, Lepto, is a general and flexible Monte Carlo (MC) to simulate 
complete lepton-nucleon scattering events and integrate cross sections. It is based on 
the leading order electroweak cross sections for the underlying parton level scattering 
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processes, and includes QCD corrections using exact first order matrix elements and higher 
orders in the leading logQ 2 parton cascade approach. The fragmentation of produced 
partons into observable hadrons is performed with the Lund string hadronization model 
0. An arbitrary configuration of a lepton and a nucleon can be defined with constraints on 
the scattering kinematics and the generated events can be transformed to different frames. 
The present version is the latest development in a series of earlier versions j|, |5| that have 
developed in stages and profitted greatly by feed-back from extensive comparisons with 
experimental data on the hadronic final state in leptoproduction. This has given valuable 
insights concerning the QCD processes of gluon radiation and boson-gluon fusion as well 
as the confinement induced hadronization of colour charged partons. The generally good 
agreement between data, in particular from the European Muon Collaboration and the 
neutrino experiments WA21 and WA25 at CERN, and the Monte Carlo shows the validity 
of the models and procedures used in the program. This can be illustrated by the results on 
longitudinal momentum spectra of different identified hadrons |J, transverse momentum 
properties, energy flows and jet structure J7j. The much higher energies available in ep 
collisions at HERA provides both qualitative and quantitative new information on QCD 
effects, such that detailed tests of MC models are now being performed [§. . 

In the following a comprehensive description is given of the theoretical framework built 
into the program (section 2) as well as the various program components (section 3) and 
their usage (section 4). 

2 Physics and MC implementation 

2.1 Kinematics 

The main kinematic relations given here are for the case of electron-proton scattering, 
e + p — > £ + H where i is the scattered lepton and H the final hadron system, but are 
of course equally valid for any lepton beam or a neutron target. Let p e , pe be the four- 
vectors of the incoming and scattered lepton, respectively, and P that of the incoming 
proton. Some basic kinematic relations are then (cf. |10[ |TT|| ) 



s = (p e + P) 2 ~4E e E p 

W 2 = ( q + py = Q* l —^ +m i 



(1) 
(2) 



x 



Q 2 = -q 2 = -(p e -pe) 2 ~4£ e £,sin 2 ^ 



(3) 




(4) 



(5) 



IP ■ q 2m v v E p (E e - E e cos 2 f ) 
P ■ q v E e -E e cos 2 | 



(6) 



y = = — 

P ' Pe V max E e 



where the '~' sign means that the lepton and nucleon masses are neglected, an approx- 
imation which is usually acceptable already at fixed target energies. Nevertheless, exact 
formulae are used in the program to avoid unnecessary approximations. Here, y/s is the 
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total invariant mass and W the invariant mass of the hadronic system H. The exchanged 
vector boson, 7/Z for neutral current interactions and W for charged current interac- 
tions, carries the momentum transfer Q 2 and the variable v is the energy of this current in 
the target rest frame. Finally, Bjorken-x and y are the convenient dimensionless scaling 
variables in the range [0, 1]. 

For the overall event kinematics there are only two independent variables and hence by 
measuring, e.g., the energy (Eg) and angle (9g) of the scattered lepton (here defined with 
respect to the incoming electron direction), all other global variables can be calculated. 
For theoretical purposes, cross section formulae etc, the variables (x, y) or (x, Q 2 ) are 
usually used as the independent ones. No assumption has here been made about the 
structure of the proton, or about the final hadronic state. By assuming the quark-parton 
model, where the current couples to a quark with four- vector = £(E p ,0,0, E p ) and 
assuming the initial and final quark to be massless one finds £ = x. The Bjorken-x 
variable can then be interpreted as the momentum fraction of the proton which is carried 
by the struck quark. However, taking QCD corrections into account this relation no longer 
holds O. 



For ep colliders like HERA where the proton beam energy is much larger than the 
electron beam energy, the phase space is somewhat special and quite elongated in the 



proton beam direction (see Fig. 2 in [10]). Although the events are certainly not evenly 
distributed in the available phase space, they are in general very asymmetric, with most 
of the final state hadrons in the 'forward' direction along the incoming proton. 

For most analyses the event kinematics need to be reconstructed. The differential 
cross sections, which are the basis of most analyses, are usually needed as functions 
of, e.g., x and Q 2 . The event kinematics can be straightforwardly obtained from the 
formulae above if the scattered lepton can be well measured. However, this is not always 
the case due to instrumental effects or because it is a neutrino. These problems are 
accentuated at ep colliders where the acceptance is limited by the beam pipes and where 
charged current interactions give an undetectable neutrino as scattered lepton. This 
implies the importance of being able to use the hadronic part of the event to reconstruct 
the kinematical variables. In the naive quark-parton model (QPM) the scattered quark 
gives the necessary information, but in reality a number of smearing effects enter (QCD 
processes, fragmentation, mass effects, jet reconstruction). Considering the hadronic final 
state as a single system, whose internal structure is of no importance, and applying energy 
momentum conservation between this system and the scattered lepton leads to methods 
such as 'Jacquet-Blondel' |L3|, [14]], 'double angle' |15j and '£' RIB| . These are particularly 



suited to ep colliders, since they minimize the effects of particles lost in the beam pipe. 



2.2 Electroweak cross sections 



In leading order electroweak theory || the differential neutral current (NC) cross section 
for the scattering of a charged lepton is (neglecting masses) given by 



d 2 a NC (e T ) Ana 1 
dx dQ 2 = ~x~Q i 



y 2 xFx(x, Q 2 ) + (1 - y)F 2 (x, Q 2 ) ± iy - y -\ xF 3 (x, Q 2 ) 



(7) 



in terms of the nucleon structure functions Fi^F^F^. Two of these are related through 
the Callan-Gross relation, 2xFi = F 2 , which holds for spin 1/2 quarks when neglecting 
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masses, intrinsic transverse momenta and order a s QCD effects. These effects are usually 
negligibly small, except at small-x [17|, [Tj|, but can optionally be included (cf. section 
2.3). More important for the structure of the electroweak theory is to take the lepton 
beam polarization into account. For a left- (L) and right-handed (R) electron one has the 
differential cross section 



d 2 cr N c(e L R p) _ 2vra 2 



(l + (1 - y) 2 ) F 2 L ' R (x, Q 2 ) + (l - (1 - y) 2 ) xF 3 L ' R (x, Q 2 )] (8) 



dx dQ 2 xQ* 
The structure functions are in the standard parton model given by 



xF 3 L ' R ( 



x,q 



53 { X( lf( x > Q 2 ) + xq/(x, Q 2 ) 
f 

53 [xqf(x, Q 2 ) - xq f (x, Q 2 ) 
f 



A L f ' R (Q 2 ) 



B L f ' R {Q 2 ) 



(9) 
(10) 



where the sum is over all flavours / and denote the quark (antiquark) density 

distributions in the nucleon (cf. section 2.4). The coefficients are given by 



.4 
B 



l Aq 2 

L f R {Q 2 ' 



e) - 2e f (v e ± a e )v f P z + (v e ± a e ) 2 (vj + a))P 2 z 
=F 2e f (v e ± a e )a f P z ± 2(y e ± a e ) 2 v f a f Pl 



(11) 
(12) 



where e/ is the electric charge (e e = —1), Vf = (T 3 / — 2e/ sin 2 Qw)l sin 29w and a/ = 
Tsf/ sin2^vK are the NC vector and axial vector couplings expressed in terms of the third 
component of the weak isospin (T 3e = — |) and the Weinberg angle 9w- Pz is the ratio of 
the Z and 7 propagators Pz = Q 2 / (Q 2 + Mg). The corresponding cross sections for left- 
and right-handed positrons, e~j^ R , are obtained from the above electron formulae by the 
replacements 

F^ R -> F R ' L , xF^ R -> -xF R ' L (13) 

The cross section for an arbitrarily polarized electron/positron beam is simply obtained 
as a linear combination of these pure left- and right-handed cross sections (cf. [|P7| ). 

The pure 7 exchange term, i.e. the one without a Pz dependence in eq. (|TID, dominates 
completely at low Q 2 , and the cross section then takes the familiar form measured in fixed 
target electron and muon beam experiments 



<i 2 cr 7 (ep) 27ra 2 



l + (l-l/) 2 )F 2 em ( 



dx dQ 2 xQ 4 
where the electromagnetic structure function is given by 

(*,q 2 )-^ 2 r- ^ -2 



x,Q< 



jpem 1 
b 2 



53 e 2 [xq f (x, Q 2 ) + xq f (x, Q 2 
f 



(14) 



(15) 



With increasing Q 2 first the 7/Z interference term (linear in Pz) and then the pure weak 
term (quadratic in Pz) become important and finally dominate the cross section at large 
Q 2 , see e.g. Fig. 3 in @. 

The differential cross sections for charged current (CC) ep interactions are given by 
d 2 acc{e~p) (1 — A)7ra 2 



dx dQ 2 

d 2 a C c(e + p) 
dx dQ 2 



Asm 4 6 w (Q 2 

{1 + X)7ca 2 
Asm 4 6 w (Q 2 



M 2 W 



M 2 W ) 



7$ 

1 1,3 

i,3 



\V Uidj Ui + (l-y) V uc 



V Utdj Ui+(l-y) V udi 



2 - 
di 



(l; 



(16) 
(17) 
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where V Ui d are elements of the Kobayashi-Maskawa matrix, and dj denote the parton 
density functions for the up-type and down-type quark flavours, respectively, and i, j are 
family indices. The beam polarization is denoted by A (±1 for a right /left-handed 
state). Considering only four massless quark flavours (u, d, s, c) and using the unitarity 

2 2 

relation J2j Vu t dj = J2j Vujd t = 1 one obtains for any lepton with fixed helicity 



d 2 g C c(£p) 
dx dQ 2 




y) 2 (d+s) for t = e L , v 
y) 2 (d + s) for £ = e£, u 
for t = e~ R , e\ 



Here, Gp = 7ra/(\/2sm 2 OwM^) is the Fermi coupling constant, Mw the H^-boson mass 
and u denotes the w-quark density u(x, Q 2 ) etc. 

The above formulae apply equally well to muon scattering, whereas for neutrino scat- 
tering some rearrangements are required. For charged current scattering, eq. ([16]) apply 
for a v with A = — 1 and eq. (|T7|) for a v with A = +1, giving the results in eq. ([18]). 
For neutral current neutrino scattering a similar correspondence applies, but the electron 
electroweak couplings must also be changed to the neutrino ones. Hence, only the pure 
Z° contribution remains in eqs. (|TTJ,[r2|) with v — Z° vector and axial vector couplings, v u 



and a v 19 



When simulating the kinematic variables according to the above differential cross 
sections one may choose the two independent variables depending on the process. The 
choice x, Q 2 for NC and x, y for CC is suitable and adopted in the program, but this can 
be changed (see LST(l) in common LEPTOU) to other combinations that may improve 
efnency, e.g., by better reflecting strong kinematic cuts. The two variables are first chosen 
according to a function that can be directly generated without any rejections, i.e. its 
primitive function can be obtained and inverted analytically. By chosing this function 
to represent the strongest variation of the cross section formulae, the remaining part can 
be taken into account by a simple rejection technique whose efficiency will be higher the 
smaller the variation in this remaining function. Formally, a function 

h v (v) = a -\ h -it + -7 (19) 

V V z V 6 

is introduced for each variable v — x,y, Q 2 , W 2 and the cross section is rewritten in the 
form (using variables x, Q 2 as a definite example) 

f da/dxdQ 2 



da = {h x (x)dx h Q2 (Q 2 )dQ 2 } a °' h a ^ (20) 



where the random variables can be chosen exactly according to the expression in the first 
bracket and the remaining factor is used for the weighting procedure. The latter simply 
means that the chosen x, Q 2 point is rejected if it results in a value for the second factor 
which is smaller than a random number (uniformly chosen between zero and unity) times 
the maximum of that factor with respect to the two variables (x, Q 2 in this case). Actually, 
any estimate larger than the maximum will do, but efficiency improves the closer to the 
true maximum it is. This estimate is, however, a constant for any fixed interaction and 
kinematic region and can therefore be obtained in the initialization phase of the Monte 
Carlo program, see subroutine LINIT. Thus new points are generated using the first factor 
until accepted by weighting with the second factor and hence the resulting events, or phase 
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space points, will have no weight associated with them (or more correctly, they all have 
unit weight). 

The parameters a, in the functions h v are in the program represented by OPTu(i) in 
common LOPTIM and are given process dependent default values in subroutine LINIT 
to optimize program speed under normal conditions. This means that the functions h v 
reflect the dominant variation of the matrix element from propagators etc. Since this 
may change with kinematics, e.g. dominance of pure Z° exchange at very large Q 2 , more 
optimal values may be found under certain conditions. 

In the initialization phase (subroutine LINIT) the simulation variables are defined and 
their effective limits calculated from applied cuts. The optimization parameters are set 



and the maximum needed for the weighting is found using an adaptation of Minuit ||20|| . 
Further, the total cross section can be obtained at this stage by numerical integration 
over the kinematic variables. In the simulation phase, see subroutine LEPTO, phase 
space points are chosen from the cross section (in subroutine LEPTOX) as discussed 
above. This Monte Carlo sampling is also used to provide an estimate of the cross section 
for the process being simulated, see PARL(24) in common LEPTOU. Since the result is 
updated with each generated event the accuracy depends on the generated statistics as 



2.3 Longitudinal structure function 



Target mass effects and the longitudinal structure function, defined by 
F L (x,Q 2 ) 



1 + I F 2 (x, Q 2 ) - 2xF 1 (x, Q 2 ) 



Q 2 



(21) 



can be included as an option (see LST(ll)) using the formalism for photon exchange. The 
cross section in eq. ( |14|) is then modified to 



d 2 a^(ep) 
dx dQ 2 



Ana 2 
~x~Q* 
2na 2 
~xQ* 



(l-V- V 2 ^-) Mx, Q 2 ) + y 2 xF 1 (x, Q 2 ) 
l + (l-y) 2 )F 2 (x,Q 2 )-y 2 F L (x,Q 2 ) 



(22) 
(23) 



where all mass effects have been absorbed into Fi which consists of the three terms 

F L (x, Q 2 ) = F? CD (x, Q 2 ) + F™(x, Q 2 ) + Q 2 ) (24) 

The QCD contribution (which is leading twist) is to order a s given by 

6n Jx y 6 7r j 1 Jx y A \ 

(25) 

The first term originates from the gluon radiation diagram and the second one from the 
photon-gluon fusion process, where the sum runs over the quark flavours (taken as u, d, s, c 
in the program). F® CD gives an important contribution at small-x, where the gluon term 



l--)yg(y,Q' 

y 
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dominates and Fl may hence be used to extract the gluon distribution fI5 |. The target 
mass correction to 0(m 2 /Q 2 ) is given by |21, |22 



Q 2 ) = 4^fx 3 f ^F 2 (y, Q 2 ) - 2^x 2 F 2 (x, Q 2 ) 



(26) 



where also the mass term in eq. ([22]) has been included. Finally the (dynamical) higher 
twist contribution to 0(1/ Q 2 ) can be written as Wl 



' l (x,Q 2 )=8—F 2 (x,Q 2 



(27) 



where the scale of this twist-4 contribution is given by the parameter k 2 (cf. PARL(19)) 
with a value around 0.03 GeV 2 obtained [22|| from SLAC data. 



In the numerical implementation of these equations the evaluation of the integrals in 
eqs. (^,^) can be performed for each event at its proper x and Q 2 -value. Since this 
tends to be time-consuming another option (see LST(ll)) is to initially (in LINIT) set 
up a grid in x, Q 2 with their values and then perform a linear interpolation to the desired 
x, Q 2 when simulating events. The latter method is preferrable and gives almost the same 
precision under normal conditions. The inclusion of Fl only affects the distribution in 
x, Q 2 , i.e. eq. (|22"D and thereby the cross section, but not the generation of the hadronic 
part of the event. For the first order QCD matrix element corrections, which are still 
generated as usual (cf. section 2.5), this is not fully consistent but should be an adequate 
approximation for inclusive properties of the hadronic final state. 



2.4 Parton density distributions 

To define the parton content of the proton for the cross section formulae above the parton 
density functions g/(x, Q 2 ), qf(x, Q 2 ) and g(x, Q 2 ) are needed. These give the probability 
to find a quark or antiquark of flavour /, or a gluon, respectively, carrying a fraction x of 
the proton momentum when probing the proton with a momentum transfer Q 2 . Several 
parametrizations of these distributions have been obtained using data, in particular from 
lepton scattering experiments, and with Q 2 - dependence according to the perturbative 
QCD evolution equations [p3fl . The fit to the data provides the x-dependence and the 
QCD parameter A. The choice among many available parton density parametrisations in 
Pythia 5.7 |24| and in Pdflib p5) is made through the switches LST(15) and LST(16) 



in common LEPTOU. 

The possibility to scatter on intrinsic charm and bottom quarks in the nucleon is 
included by an option, see LST(15). The hypothesis of intrinsic charm quarks in the 
proton was introduced |26| as an attempt to understand a large discrepancy between early 
charm hadroproduction data and leading order perturbative QCD (pQCD) calculations. 
This discrepance has largely disappeared as more data have been collected and next-to- 
leading order (NLO) pQCD calculations have been made. Still, however, there are certain 
aspects of charm production data which are difficult to understand within the pQCD 
framework, but are natural if the intrinsic charm hypothesis is basically correct (see |27| 



and references therein) . Intrinsic charm (IC) corresponds to a Fock-state decomposition of 
the proton wave function, \p) = a\uud)+j3\uudcc) + ..., with a small, but finite, probability 
j3 2 (PARL(12)) for an intrinsic cc pair as a quantum fluctuation. The normalization of 
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this component is the key unknown, although it should decrease as I/tuq. Originally, a 
1% probability was assumed, but later investigations based on EMC data on the charm 
structure function F£(x, Q 2 ) [28| indicate a somewhat smaller but non-vanishing level; 



0.3% H and (0.86 ± 0.60)% |30|. 

From the IC model one obtaines an effective charm quark density p7 | 



c w (x) = /5 2 1800a; 2 



-(1 - x)(l + 10x + x 2 ) + 2x(l + x) \nx | 



which gives a charactersitic hard momentum spectrum with (x c ) = 2/7. The Q 2 depen- 
dence from normal leading log GLAP equations have been calculated for IC p^ , but 
can be taken into account through a simple extension of the parameterisation in Eq. (|28|) 
PTfl . This quark density is then used in the electroweak cross section fomulae to simu- 
late the scattering on such an intrinsic charm (or bottom) quark, with parton showers 
and hadronization added as usual. The scattered (anti)charm quark gives a charmed 
hadron in the current fragmentation region and the remaining partner (anti) charm quark 
in the proton remnant gives a charmed hadron in the target fragmentation region. For 



phenomenological studies with this model, see p7|. 



2.5 First order QCD processes 

The leading order parton level process is V*q — > q, where V* is the exchanged virtual 
boson j/Z or W. In first order QCD the gluon radiation or QCD Compton process, 
V*q —>■ qg, and the boson-gluon fusion (BGF) process, V*g — ► qq, appear and can be 
included with their matrix elements (ME) |3"T| , [3^, [33], |34[ . Quark masses are not included 
in these ME, but a threshold factor is applied for boson-gluon fusion into heavy quarks. 
For more accurate heavy quark simulations the AROMA Monte Carlo [3_5[ may be used. 



The first order ME are rather complicated functions involving three new degrees of 
freedom corresponding to energy, polar angle and azimuthal angle of one final parton 
(the other is then determined by energy-momentum conservation). In terms of the more 
suitable variables [3~2"| x p = x/l; and z q = P ■ p q /P ■ q, where £ is the momentum fraction 



that the incoming parton take of the proton and p q the four-momentum of, e.g., the final 
quark, the cross sections are five-fold differential 

(29) 



dx dQ 2 dx p dz, 



Here is the parton azimuthal angle with respect to the lepton scattering plane (in the 
7P cms) and its distribution has the form 

da = d<TQ + dai cos <fi + da 2 cos 20 (30) 

where dai depends on the other four variables. When averaging over 0, only the first 
term contributes and the 0-dependence can therefore often be neglected, but in dedicated 
analyses it can be observed |36[ . The invariant mass squared of the two emerging partons 



(qg or qq) is labelled s. 

The matrix elements have soft and collinear divergences that may be partly cancelled 
by virtual corrections and partly absorbed in the parton density functions. The QCD- 
Compton cross section a qg diverges as 1/(1 — x p )(l — z q ) and the boson gluon fusion a qq 
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as l/z q (l — z q ). To avoid these singularities in a Monte Carlo simulation procedure, it is 
common practice to impose a cut-off on the matrix elements. In Lepto, several different 
cut-off procedures are available through LST(20). In the W scheme it is required that 
each pair of emerging partons (including the proton remnant) have a minimum invariant 
mass, expressed as = (pi + Pj) 2 > y cu tW 2 . As a new default in the present version 



of the program, we use a variation of the 'mixed scheme' ||34|| . In this l zs > scheme the 
applied cuts are z q , min < z q < l—z^ min and s > s min , with z q , min and s min set by PARL(8) 
and PARL(9), respectively. As discussed in |j37| , this gives some advantages compared to 



previously used cut-off schemes. In particular, it gives a better division of the available 
phase space into a region for hard emission, which is best described by exact matrix 
elements, and the soft and collinear regions where higher orders are important and can 
be taken into account by leading log parton showers (cf. section 2.6). 

To decide on an event-by-event basis whether to generate one of the first order event 
types (qg- or gg-event) rather than the leading order process (g-event), the probability for 
each event type must be available as a function of the kinematic variables x, Q 2 . The q- 
event probability is taken as P q = 1 — P qg —P qq and the probabilities P qg and P qq are defined 
as the integral of the relevant first order matrix element over the above three variables, 
divided by the overall differential cross section da/dxdQ 2 . The integration over is trivial 
and the variable z q can be integrated analytically leaving only the x p -integration, which 
involves parton density parametrizations, to be performed numerically using an adaptive 
Gaussian method. If the cut-off is chosen too low the calculated probabilities for the first 
order processeses P qg + P q q can be larger than unity at the given x and Q 2 . In this case the 
cutoff is increased, the y cut in the W scheme and the z min or the s min in the zs scheme, 
until 1— PARL(13)< P qg + P qq < 1. The effective cut-off actually used in a generated 
event is stored in PARL(27). 

To save computer time the probabilities for qg- and gg-events can be calculated at the 
initialization time (LINIT) and stored on a 'grid' in the x, W or x, y plane. In the following 
event generation phase the probabilities at any x, Q 2 (or x, y) point, chosen according to 
section 2.2, are then obtained by linear interpolation on this grid. The switch LST(19) 
regulates the use of such grids and provides an option with a grid that automatically 
adjusts to the kinematical region chosen by the user. There is also a possibility to calculate 
these integrals for each event, and thereby avoid interpolation errors. 

Having chosen to generate a first order event based on the probabilities P q , P qg and 
P qq , the internal variables x p , z q and <ft are generated in turn from the matrix element 
formulae where the subsequent variables are integrated out. Given the values of all five 
variables the four-momenta of the scattered lepton and partons can be calculated. 



2.6 QCD parton shower evolution 

In order to take higher than first order QCD effects into account the parton shower 
(PS) approach has been implemented as described in detail in ref. [|l2|]. This has the 
advantage that arbitrarily high orders in a s can be simulated, but only in the leading 
log Q 2 approximation as opposed to the exact treatment in fixed order ME. Higher order 
effects are important at high energies where multiple parton emission can give rise to 
multijet events as well as affect the internal properties, such as hardness and width, of a 
jet [T3|, |5S| and the overall structure of the event in terms of, e.g., energy flows. 
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In DIS the quark struck by the electroweak boson can emit partons both before and 
after the boson vertex giving rise to initial and final state parton showers, respectively. A 
parton close to mass-shell in the incoming nucleon can initiate a parton emission cascade 
(or shower) where in each branching one parton becomes increasingly off-shell with a 
space-like virtuality (m 2 < 0) and the other is on-shell or has a time-like virtuality (m 2 > 
0). This initial state space-like shower results in a space- like quark which interacts with 
the electroweak boson that turns it into an outgoing quark which is either on-shell or 
has a time-like virtuality In the latter case a final state, time-like shower will result 
where the off-shell mass is reduced by branching into daughter partons with decreasing 
off-shell masses and decreasing opening angles. This shower continues until all partons 
are (essentially) on-shell. Any parton with a time-like virtuality from the initial state 
shower will develop similarly. The general behaviour of initial and final state showers 
are similar since they are both based on the branching processes q — > qg, g —>■ gg and 
g —>■ qq as described by the GLAP equations [23 in the leading logQ 2 approximation of 
perturbative QCD. 

The final state radiation is analogous to parton radiation in e + e~ — > qq and is theo- 
retically well developed and tested against data. The routine LUSHOW in Jetset |2j[] 



is therefore used for all time-like showers. The evolution is based on the Sudakov form 
factor, which expresses the probability that a parton does not branch between some initial 
maximum virtuality and some minimum value. From this one can find the mass of the 
decaying parton, the energy fractions in the branching and the flavours of the daughter 
partons. The process is iterated with a reduced virtuality until all parton virtualities are 
below some cutoff m 2 , around 1 GeV 2 . The technical details are given in |12| , but it should 



be noted that coherence in soft gluon emission is taken into account through angular or- 
dering (decreasing opening angles in subsequent branches) and that p 2 ± ~ z(l — z)m 2 is 
used as argument in a s as suggested by studies of coherence effects. 



The initial state radiation is performed using the 'backwards' evolution scheme f39 
where the shower is constructed from the hard electroweak interaction backwards with 
decreasing virtualities down to the on-shell parton from the incoming nucleon. A modified 



version JT2[ of the routine PYSSPA in Pythia |^0[ is used for the space-like shower. This 
is a more complicated process, e.g. since the nucleon parton density functions must be 
taken into account (which tend to reduce the amount of radiation). In addition it is not 
so well tested by the more 'messy' hadron collision data. When combining the initial 
and final state radiation to get the complete model, special precautions have been taken 
T2| to preserve energy-momentum conservation and keep the normal definitions of the 
kinematic variables for the electroweak scattering such that they will be obtained from 
the scattered lepton as usual, in particular that Bjorken-rr is preserved. 

The amount and hardness of the initial and final radiation depends on the off-shellness 
of the struck parton just before and the partons after the boson vertex. These virtualities 
are chosen, using the Sudakov form factor, between the lower cut-off and a maximum 
value to be given by the energy or momentum transfer scale in the process. The default 
procedure is to add the parton shower to a matrix element event. The possible final states 
are then q, qg and qq events. In the case of a q event the maximum scales are set by the 
matrix element cut-off (e.g. i/cutW 2 in the M^-scheme) since emissions harder than the 
cut-off would be double-counting. In the case of a qg or qq event the maximum virtuality 
scale for the final state shower is set to s. As maximum virtuality for the backwards initial 
state shower it is natural to use the mass-squared of the quark propagator just before the 
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boson vertex. This is given by the known four-vectors of the exchanged boson and the two 
final partons from the matrix element, but depending on the underlying Feynman diagram 
in the amplitude different combinations are possible leaving some remaining ambiguity. 
The largest of these possible virtualities is being used (cf. subroutine LSCALE). 

In the case where the parton shower is used without the matrix element there is a large 
ambiguity in which maximum scale should be used. Whereas only one scale is present in 
e + e~, any function of Q 2 and W 2 may be possible in DIS. Although Q 2 and W 2 can often 
be of similar magnitude, cf. eq. (2), at small x- values W 2 is much larger than Q 2 giving 



rise to drastically different amounts of radiation |41| . Given this uncertainty, different 
scales are available (cf. LST(9)) based on different motivations. The phase space limit 
is given by W 2 , such that a smaller scale (Q 2 ) would cut off the tail of high-p^ parton 
emission. Refering to the first order ME parton level one finds ~ Q 2 (l — x) for x — ► 1 
and (p 2 ± ) ~ Q 2 ln(l/x) for x — > |3l|. A suitable interpolation between these limiting 
behaviours is the choice Q 2 {1 — x)max(l, In -) which is the default choice when the matrix 
elements are not used. 

The parton shower approach has some shortcomings due to its approximate nature. 
The separation of initial and final state parton emission implies the neglect of interference 
terms between the two and is not gauge invariant. The use of the leading logarithm 
approximation means that the emission of partons that are soft or close to the directions 
of the emitting partons should be well described, while the emission of hard partons at 
large angles could be mistreated. Therefore, the rate of events with extra hard partons 
that give rise to separate jets, i.e. multijet events, need not be well accounted for. The 
use of matrix elements is preferable for these purposes and is therefore the default option. 
Starting with the ME the hard emission is generated and extra, but softer emissions are 
then added using PS. An advantage of this procedure is, of course, that the hard parton 
emission is properly treated using the first order ME. 



2.7 Nucleon remnant and hadronization 

The remnant system is the target nucleon 'minus' the parton entering the hard scattering 
system (initial parton showers and matrix elements). This interacting parton can be either 
a valence quark, a sea-quark or a gluon. 

When the interacting parton is a valence quark the nucleon remnant is simply a diquark 
composed of the two left-over valence quarks as spectators. In the Lund model || a colour 
triplet string is stretched between the colour triplet charged struck quark and the diquark 
which is a colour antitriplet. This system is then hadronized in the usual way H |23J by 
the production of quark-antiquark and diquark-antidiquark pairs from the energy in the 
field, leading to hadron production. The proton remnant diquark is not a single entity; 
its two quarks may go into a leading baryon but they can also be separated to produce a 
leading meson followed by a baryon. 

In case the interacting parton quark (q s ) or antiquark the nucleon remnant 

contains the corresponding antiquark or quark in addition to the three valence quarks 
(q v ). This more complicated four-quark system q v q v q v q s or q v q v q v q s must be taken into 
account to conserve the flavour quantum numbers. 

In the conventional way (default in Lepto version 6.2 and earlier) the following treat- 
ment has been used. If q s = u or d it is cancelled against a corresponding valence quark 
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leaving a simple diquark system to be treated as above. For other flavours of q s it is joined 
with a valence quark of arbitrary flavour into a meson (M = q v q s )- The q s is assumed to 
have no specific dynamic properties such that this splitting process into a meson and a 
diquark should be similar to normal hadronization. The meson is then given a fraction 
z of the remnants energy-momentum (E + p z ) along the beam direction from a probabil- 
ity distribution P{z) (cf. LST(14)) and only a small Gaussian p± (cf. PARL(14)). The 
left-over diquark, with longitudinal momentum given by 1 — z and equal but opposite p±, 
forms a string system with the scattered quark and hadronization proceeds as usual. If an 
antiquark (q s ) was scattered the remnant is a four-quark system q v q v q v q s which is treated 
similarly to the previous case. Here, the corresponding quark (q s ) is combined with a 
random diquark giving a baryon (B = q v q v q s ) leaving the remaining valence quark to 
form a string system with the scattered antiquark. The split of the remnant is as before, 
taking account of the masses in the distribution for z (cf. LST(14)). 

In Lepto 6.3 a modified treatment of sea quarks in the remnant was introduced which 
is now default (cf. LST(35)). The essential difference is that the sea quark partners (q s ) 
are treated dynamically and also u and d quarks can be considered as sea quarks. The 
interacting quark is assigned to be a valence or sea quark from the relative size of the 
corresponding parton distributions q v (%i,Qi) and q s (xi,Ql), where X\ is the momentum 
fraction of the quark 'leaving' the proton and Q\ is the relevant scale (typically the cutoff 
Qq of the initial state parton shower). In case of a valence quark the previous treatment 
is used, but in case of a sea quark a new treatment is used. The left-over partner q s 
is given a longitudinal momentum fraction from the Altarelli-Parisi splitting function 
P{9 ~~ > q<i) an d the transverse momentum follows from the masses of the partons in 
the splitting. Essentially the same results are obtained if the longitudinal momentum 
fraction is chosen from the corresponding sea quark momentum distribution. The former 
approach is presently used since this allows the mechanism to be simply implemented in 
the initial state parton shower routine as an additional, but non-perturbative, g — > qq 
process. This partner sea quark will then be at the end-point of a string and not, as 
previously, go directly into a hadron together with another spectator parton. Depending 
on the momentum of the partner sea quark, this new string may extend more or less into 
the central region and through hadronization contribute to the particle and energy flow 
in the forward region. In particular, the transverse forward energy flow will be enhanced 



42, 37 and improve the agreement with HERA data. 



In boson-gluon fusion the removed gluon leaves the three valence quarks in a colour 
octet state. This remnant is split into a quark and a diquark, chosen with random flavours, 
which form two separate strings with the antiquark and quark, respectively, produced in 
the fusion process. Again the split of the remnant involves the same longitudinal mo- 
mentum sharing and a Gaussian transverse momentum. For the order a s gluon radiation 
process (gg-event) the string is stretched from the scattered quark via the gluon to the 
target remnant. 

In the parton shower case, the backwards evolution always results in one parton being 
removed from the nucleon as in the above cases such that the same procedures can be 
applied. The additional partons emitted in the PS case will, however, lead to a more 
complicated string configuration. The string follows the colour flow of the parton shower 
such that it starts from a colour triplet quark and goes via a number of colour octet 
gluons, which are kinks on the string, before ending up on a colour antitriplet antiquark 
or diquark. Where quark-antiquark pairs have been produced in the shower, the colour 
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flow will be broken resulting in a termination of the first string piece and the start of 
a new one. The string system may thus be divided into subunits which then hadronize 
separately. 

The ME and PS emissions may give a varying number of soft or collinear partons, 
depending on the details of the cut-offs. Although such partons cannot be observed as 
separate jets, they may give a 'softening' and 'fattening' of jets. The Lund string model 
is particularly suitable in this context, since it provides a stability in the sense that the 
hadron level result will not depend strongly on the presence of extra soft partons. Rather, 



one obtains a smooth transition to a configuration without them p3, 44|. The independent 



hadronization model, available as an option in [24] does not have the same property and 
is therefore not recommendable. 

In this context one should also note that the two-string configuration for sea-quark 
initiated processes provides a desirable continuity between the two-string gluon-initiated 
gg-events and the one-string quark-initiated g-events. Depending on the partner sea-quark 
momentum, the corresponding string will extend more or less into the central region in 
rapidity. The hadronization of this extra string will contribute to the particle multiplicity 



and energy flow in this region [42, 37 



The parameters for the hadronization process in Jetset [^J] are obtained from fits to 
e + e~ data and are assumed to be the same in DIS based on jet universality. Nevertheless, 
they depend on which QCD effects are explicitly included in the Monte Carlo simulation. 
The default values are suitable when higher orders are taken into account via parton 
showers, whereas with first order ME alone the hadronization should be made slightly 
'softer' and 'wider' to account for the additional parton emission not simulated explicitly. 



2.8 Soft colour interactions and rapidity gaps 



The rapidity gap events discovered in deep inelastic scattering at HERA [45 are usually 



interpreted in terms of pomeron exchange models [J46|| . Although this seems to work 
reasonably well phenomenologically, there is no satisfactory understanding of the pomeron 
and its interactions mechanism. As an alternative, we have introduced a model E2L 



based on soft colour interactions (SCI) that give rise to rapidity gap events without using 
the concept of a pomeron. 

At small Bjorken-x (10 -4 — 10~ 2 ), where the rapidity gap events occur, the events 
are frequently initiated by a gluon from the proton. This can either be directly from the 
boson gluon fusion matrix element or after the initial state parton shower, including a 
possible split in the sea quark treatment. In the conventional string hadronization model 
this gives two separate strings from the q and q to the proton remnant spectator partons, 
where the gluons from the parton shower are kinks on the string, thereby causing particle 
production over the whole rapidity region in between. The new hypothesis introduced 
here is that additional non-perturbative soft colour interactions may occur. These have 
small momentum transfers, below the scale Ql defining the limit of pQCD, and do not 
significantly change momenta from the perturbative phase. However, SCI will change the 
colour of the partons involved and thereby change the colour topology as represented by 
the strings. Thus, it is proposed |42j that the perturbatively produced quarks and gluons 



can interact softly with the colour medium of the proton as they propagate through it. 
This should be a natural part of the processes in which 'bare' perturbative partons are 
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'dressed' into non-pertubative quarks and gluons and the formation of the confining colour 
flux tube in between them. 

Lacking a proper understanding of such non-perturbative QCD processes, a simple 
model is used to describe and simulate these interactions. All partons from the hard 
interaction (electroweak + pQCD) plus the remaining quarks in the proton remnant con- 
stitute a set of colour charges. Each pair of charges can make a soft interaction changing 
only the colour and not the momenta, which may be viewed as soft non-perturbative 
gluon exchange. As the process is non-perturbative the exchange probability cannot be 
calculated so instead it is described by a parameter, PARL(7). The number of soft ex- 
changes will vary event-by-event and change the colour topology of the events such that, 
in some cases, colour singlet subsystems arise separated in rapidity. In the Lund model 
this corresponds to a modified string stretching and rapidity gaps may arise when a gap 



at the parton level is not spanned by a string, as illustrated in [42]. In particular, when 
the hard process starts with a gluon from the proton, leaving a colour octet remnant and 
giving a colour octet hard scattering system, a soft colour exchange between the two octet 
systems can give two colour singlets separated in rapidity. 

SCI may thereby give a string system, including a valence diquark from the proton 
remnant, which has a small invariant mass. Such systems are not optimally treated 
in Jetset regarding production of one- or two-particle final states and taking isospin 
constraints into account. We have therefore constructed a new treatment (implemented 
in subroutine LSMALL) for such systems giving two particles if kinematically possible and 
otherwise one. When one particle is made, the isospin restriction from the effective isospin 
singlet exchange in SCI is taken into account to prevent A production. The invariant mass 
of such small systems is roughly given by the constituent masses of the partons in the 
proton remnant and the transverse momentum in the remnant split. With the default 
contituent masses (m q = 325 MeV and m qq = 650 MeV) and width of the transverse 
momentum in the remnant split (350 MeV), the invariant mass of the consituents in 
the proton remnant exceeds the mass of the proton itself. Therefore the possibility has 
been introduced to reduce the constituent masses of the remnant partons when the target 
remnant is more complicated than a simple diquark (cf. PARL(20)). 

SCI also gives rise to other string topologies without gaps, but where the string goes 
back and forth in rapidity between the partons. In such cases, there will be more particles 
and energy per unit of rapidity. This contributes to a better description of, e.g., the 



forward energy flow observed at HERA as discussed in [37]. 



3 Description of program components 

The present program, Lepto version 6.5, is a backwards compatible update of versions 
6.1 H to 6.4. It is a further development of earlier versions [|j], but is not backwards 
compatible with those due to the changed conventions in Jetset version 7 now being 
used. Previous knowledge of earlier versions is not necessary, but is helpful since the main 
structure and use of the program has been kept similar. The code is written completely in 
standard FORTRAN77 and should therefore run on any computer with such a compiler. 
Single precision is normally used, but double precision is being used when required. 
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3.1 Subroutines and functions 



The following routines should or may be called by the user: 



SUBROUTINE LINIT(LFILE,LEPIN,PLZ,PPZ,INTER) 



Purpose: 

Arguments: 
LFILE : 

= 
< 
> 
Remarks 



LEPIN : 



PLZ, PPZ : 



INTER 



Procedure: 



to initialize the event generation procedure and, optionally, integrate cross 
section. 

logical file number containing weights for first order QCD, see LST(8). 
the weights are calculated but not save on file, no file is used, 
the weights are calculated and stored on file number -LFILE. 
the weights are read from file number LFILE. 

When using weights from a file, the same conditions (interaction, cuts etc) 

must hold as when the weights were calculated. The relevant quantities 

are checked and a mismatch results in an error message, cf. LST(3). There 

is no strong reason, in this upgraded program, to store weights on a file 

since their calculation is fast enough to be repeated in each run. 

type of lepton, i.e. 11 = e~, 12 = u e , 13 = 14 = and negative values 

for the coresponding antiparticles, i.e. Jetset |24j] code. 

momentum (GeV/c) for incoming lepton and nucleon, respectively, along 

the z-axis (if both non-zero, i.e. colliding beams, they must have opposite 

signs). Colliding beams which are not along a common axis, or variable 

beam energies, are possible using LST(17). 

type of interaction to be simulated. 

electromagnetic (EM), i.e. 7 exchange. 

weak charged current (CC), i.e. W ± exchange. 

weak neutral current, i.e. Z° exchange. 

neutral current (NC), i.e. 7/Z exchange. 

Various constants are calculated, effective limits on kinematic variables 
are derived from the cuts and, depending on the interaction, integration 
variables are chosen and parameters are set. For the Monte Carlo rejec- 
tion technique to give unweighted events, the maximum of the differential 
cross section (having factored out a suitable function that can be sim- 
ulated exactly, see section 2.2) is found using an adaptation of Minuit 
The cross section is calculated using numerical integration (LST(IO)) 



20 



and stored (PARL(23)). Optionally, grids with probabilities for first order 
QCD processes are set up (LST(8)) and grids for including the longitudinal 
structure function (LST(ll)). 

SUBROUTINE LEPTO 

Purpose: to administer the generation of one event of the kind specified by the last 
LINIT call. 

Procedure: Beam energies are taken from LINIT (optionally from LUJETS, see LST(17)). 

A phase space point x, Q 2 is chosen according to the differential cross 
section (using LEPTOX), and the parton level system, optionally with 
QCD corrections (cf. LST(8)) is set up. Hadronization and decays are 
performed, via LUEXEC, and the event is transformed to the selected 
frame (LST(5),LST(6)). The Monte Carlo estimate of the cross section, 
PARL(24), is updated with each event. 



16 



Lepto 6.5 manual 



Remarks: Under some conditions, an error may have occured and the event should 
be rejected, see LST(21). 

SUBROUTINE LFRAME(IFRAME,IPHI) 



Purpose: 
Arguments: 
IFRAME : 



to transform the event between different frames. 

specifies the desired frame (as for LST(5)). 
= 1: hadronic CM frame, z-axis along exchanged boson. 
=2: lepton-nucleon CM frame, z-axis along lepton. 
=3: lab system as defined by last user call to LINIT. 
=4: as 3, but z-axis along exchanged boson. 
IPHI : specifies whether to include a random rotation for the azimuthal angle, 0, 

of the lepton scattering plane. 
=0: no rotation, scattering plane is x — z plane. A possible earlier rotation is 
undone. 

= 1: random rotation in 0, made in lepton-nucleon CM frame. Not made for 
IFRAME=1. 

Remark: The present frame is stored in LST(28), LST(29) and is updated by LFRAME. 

Transforming the event with other routines can therefore cause errors in a 
following LFRAME call. 

SUBROUTINE LPRWTS(NSTEP) 

Purpose: to print a table of the QCD weights in common LGRID, i.e. the probabil- 
ities for q-, qg- and gg-events in first order QCD. Only the values on each 
NSTEP point in the x,W grid is printed, i.e. NSTEP=1 prints all grid 
points. 

SUBROUTINE LWBB(ENU) 

Purpose: to give energy (ENU in GeV) of a (anti-)neutrino chosen from a sim- 
ple parametrization (defined by DATA statement within the routine) of a 
wide band beam. The energy is actually chosen from the beam spectrum 
weighted with the energy to take into account the linear rise with energy 
of the cross section, which Lepto does not account for. This routine is a 
simple example to be replaced by a more realistic beam energy distribution 
if necessary. 

Remark: For running with variable beam energies, see LST(17). 
SUBROUTINE LTIMEX(TIME) 

Purpose: to get cpu execution time (in seconds) since start of job. Interface to 
machine dependent routine, by default TIMEX (Z007 in the CERN library) 
which can simply be changed, or TIME set to zero since timing information 
is not essential (although useful). 

SUBROUTINE LNSTRF(X,Q2,XPQ) 

Purpose: to give the parton distribution functions per nucleon for a nucleus defined 
by PARL(l) and PARL(2). Arguments as in LYSTFU which is called. 
This routine is called internally, but can also be used separately after ini- 
tialization by LINIT. 



SUBROUTINE LYSTFU(KF,X,Q2,XPQ) 



Lepto 6.5 manual 



17 



Purpose: to give the parton distribution functions for the target hadron through an 
interface to subroutine PYSTFU in Pythia 5.7 complemented with 
options for intrinsic charm and beauty quarks in the nucleon. 

particle flavour, e.g. 2212=p, 2112=n, -2212=p, -2112=n. 
momentum fraction carried by the parton. 
momentum transfer scale Q 2 - 

array (—6:6) that on return contains the values of the momentum-weighted 
parton densities, i.e. x ■ q(x,Q 2 ) and x ■ g(x,Q 2 ). Index: 0=g, l=d, 2=u, 
3=s, 4=c, 5=6, 6=t and —l=d, —2—u, — 3=s, — 4=c, —5=6, —6=i. 
Different parametrizations of structure functions are implemented and can 
be selected by LST(15) and LST(16). This routine is called internally but 
can also be used separately after initialization by UNIT. 

In the following list all subroutines (S) and functions (F) are briefly described. The 
order is as they appear in the code and reflects the flow in the program. Further details 
about their purpose and procedures used are given by comments in the code. All routine 
names start with characteristic letters to indicate origin and avoid name clashes. L, 
or D for real functions, is for Lepto routines in general; FL for routines related to 
longitudinal structure function; LY for modified PYTHIA routines (20 
Minuit routines 



Arguments: 
KF : 
X : 
Q2 : 
XPQ : 



Remarks: 



LM for modified 



20j; GADAP, RIW and DV for routines related to different integration 



procedures. 



Routine Purpose 

LTIMEX (S) to give execution time since start of job, see above. 

LEPTOD block data to give default values to all switches and parameters. 

UNIT (S) to initialize program package, see above. 

LEPTO (S) to administer the generation of an event, see above. 

LEPTOX (S) called by LEPTO to generate kinematic variables, within the applied cuts, 
according to the differential cross section giving unweighted events. Op- 
tionally include F L . Update Monte Carlo estimate of cross section, PARL(24), 
select lepton helicity and target nucleon. 

LKINEM (F) called by LEPTOX to calculate various kinematic variables and optionally 
(LST(2)) to reject event if outside given kinematic limits in CUT. 

LQCDPR (S) called by LEPTO to obtain probabilities for first order QCD processes, 
either by linear interpolation of weights stored on a grid set up in LINIT, 
or direct calculation event-by-event, see LST(19). 

LQEV (S) called from LEPTO to generate parton system for g-event, i.e. without 
QCD processes. 

LQEVAR entry in LQEV used by Ariadne [H7] for some boson-gluon fusion events. 
LQGEV (S) called from LEPTO to generate parton system for gg-event, i.e. first order 
QCD gluon radiation. 

LQQBEV (S) called from LEPTO to generate parton system for gg-event, i.e. first order 

QCD boson-gluon fusion. 
LXP (S) called from LQGEV and LQQBEV to generate a value of x p from the 

relevant QCD matrix element folded with parton density functions (for 

given x, Q 2 , but z q and integrated out). 
LZP (S) called from LQGEV and LQQBEV to generate a value of z q , from QCD 

matrix element (for given x,Q 2 ,x p , but (ft integrated out). 
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LQMCUT (F) to apply cuts on QCD parton configuration to take quark mass effects into 

account and ensure ability to apply string hadronization. 
LAZIMU (S) called from LQGEV and LQQBEV to generate the azimuthal angle <ft, of 

parton plane w.r.t. lepton scattering plane, from first order QCD matrix 

elements eq. fl30|), for given x, Q 2 , x p , z q . 
DSIGMA (F) differential cross section da/dx p for first order QCD processes. 
DSIGM2 (F) modified DSIGMA with a variable substitution to speed up integration 
DQCD (F) first order QCD differential cross sections da/dx p dz q (section 2.5) 
DQCDI (F) first order QCD differential cross section da/dx p obtained after analytical 

integration over z p and factoring out 1/(1 — x p ). 
LFLAV (S) selects flavour of struck quark and outgoing quark (with flavour mixing 

in CC) and defines the corresponding parton system from the nucleon 

remnant, also applies threshold factor for charm and heavier quarks. 
LREMH (S) gives energy-momentum fraction for target remnant split, cf. LST(14); also 

selects remnant flavours (not for PS case). 
LPRIKT (S) generates magnitude and azimuthal angle for Gaussian primordial k±_ of 

parton in nucleon, see PARL(3). 
LFRAME (S) transforms event to different frames, see above. 
LWBB (S) selects energy from a neutrino wide band beam, see above. 
LWEITS (S) integrates first order QCD matrix elements to set up grid in x, W or x, y 

of weights for gluon radiation and boson-gluon fusion and finds maxima 

used for QCD simulation, see LINIT and common LGRID. 
LPRWTS (S) prints QCD weights, see above. 

LSIGMX (S) called by Minuit routines to calculate differential electroweak cross sec- 
tion divided by the optimization function to obtain the maximum for the 
weighting procedure. 

LXSECT (S) called from LINIT to integrate cross section using Gadap, Riwiad or 

Divonne, see LINIT and LST(IO). 
RIWIBD (S) substitutes block data for Riwiad |48[| , print flag changed. 
DVNOPT (S) substitutes block data for Divonne g[J, print flag changed. 
DFUN (F) integrand for Divonne integration, calls RIWFUN. 
RIWFUN (F) integrand for Riwiad integration, calls DCROSS. 

DCROSS (F) differential electroweak cross section used as integrand for numerical inte- 
gration, see LXSECT. Integration variables given by LST(l), function is 
zero outside the cuts in array CUT. 

DLOWER, DUPPER (F) lower and upper limit on second kinematic variable (y,Q 2 or 
W 2 ) given a value of the first variable (x) and the cuts. Used together 
with DCROSS for cross section integration. 

called from LINIT to tabulate integrals in x, Q 2 -grid for the longitudinal 
structure function, see LST(ll). 

called from LEPTOX to obtain longitudinal structure function from inter- 
polation on grid from FLTABL, see LST(ll). 

called from LEPTOX to obtain longitudinal structure function by integra- 
tion event-by-event, see LST(ll). 

integrand for quark contribution to QCD longitudinal structure function, 
integrand for gluon contribution to QCD longitudinal structure function, 
integrand for target mass correction, see LST(ll). 
generates soft colour interactions (SCI), see section 2.8. 
LECSWI,LEASWI (S) switches colour and anticolour pointers for SCI. 



FLTABL 


(S) 


FLIPOL 


(S) 


FLINTG 


(S) 


FLQINT 


(F) 


FLGINT 


(F) 


FLTINT 


(F) 


LSCI 


(S) 
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LSMALL (S) 
LSHOWR (S) 



creates hadrons from small mass colour singlet systems including diquarks. 
administer parton cascade evolution added to a g-event. 
LMEPS (S) administer parton cascade evolution added to a qg- or gg-event from first 
order QCD matrix elements. 

gives scale for maximum virtuality in parton showers, see LST(9). 
simulate initial state parton cascade evolution, modification Jl2| of routine 
in M- 

LYREMN, LYSPLI (S) treatment of target remnant and primordial kj_ when using par- 



LSCALE 
LYSSPA 



(S) 
(S) 



ton showers, modifications |12| of routines in 
LMCMND,LMINTO,LMIDAT,LMINEW,LMPRIN,LMPINT,LMRAZZ,LMSIMP modified 

Minuit routines [EDI to find maximum of cross section. 
GADAP, GADAP2, GADAPF one- and two-dimensional adaptive Gaussian integration 
routines 



50]. 



LNSTRF (S) gives parton distributions in a nucleus, see above. 

LYSTFU (S) interface to subroutine PYSTFU in Pythia 5.7 |24[ to give parton distri- 
bution functions, see above. 



3.2 Common blocks 

Most of the communication between the user and the program is via the switches and 
parameters in the common blocks. The user need mainly be concerned with common 
LEPTOU since all others are essentially for internal use. All variables are given sensible 
default values in block data LEPTOD, as shown by (D=...) below. These values may 
be changed by the user to modify the behaviour of the program. Note, however, that 
this should usually be done before calling LINIT and that some of the parameters are 
interrelated. Variables whose name start with D are in double precision. The generated 
event is stored in common LUJETS, described in 



COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U 

Purpose: contains input switches (LST(1)-LST(20),LST(34),LST(35)) and input pa- 
rameters (PARL(1)-PARL(20)) to specify physics, kinematic cuts and nu- 
merical procedures, as well as output flags (LST(21)-LST(40)) and output 
variables (PARL(21)-PARL(30)). Overwriting default values should be 
made before calling LINIT. 



Parameters: 
CUT(l), CUT(2) 
CUT(3), CUT(4) 
CUT(5), CUT(6) 
CUT(7), CUT(8) 



(D= 10 , 1.) lower and upper limit of Bjorken-x variable. 
(D= 0., 1.) lower and upper limit of y variable. 
(D= 4., 10 8 ) lower and upper limit of Q 2 (GeV 2 ). 
(D= 5., 10 8 ) lower and upper limit of W 2 (GeV 2 ). 
CUT(9), CUT(IO) : (D= L, 10 s ) lower and upper limit of variable u (GeV). 
CUT(ll), CUT(12) : (D=l., 10 8 ) lower and upper limit of scattered lepton energy (GeV, 

in frame defined by LINIT call). 
CUT(13), CUT(14) : (D=0., 3.1416) lower and upper limit of lepton scattering angle (in 
radians), with respect to incoming lepton in frame defined by LINIT call. 
Remarks: These cuts are applied already when choosing kinematic variables, before 
evaluating cross section formulae and structure functions, and will there- 
fore be more efficient than applying cuts later on in the users program. The 
cross section estimates take these cuts into account. CUT(11)-CUT(14) 
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are less efficient, since they refer to a special frame and the user should, if 
possible, translate them into cuts in x,y,Q 2 . 

(D=0) choice of the two independent variables to be used for simulation 

and numerical integration of cross section. 
0: program makes a 'best' choice for simulation efficiency depending on which 

process is to be generated (1 if 7-exchange included, else 2) cf. LST(31). 
1: x and Q 2 
2: x and y 
3: x and W 2 

(D=l) choice of simulation and applying cuts in array CUT. (Used inter- 
nally also with negative values.) 
1: kinematic variables, LST(l), chosen from differential cross section and cuts 
applied. 

2: variables (x, y) supplied by user via LEPTOU, cuts applied. 
3: variables (x, y) supplied by user via LEPTOU, cuts not applied. 

(D=5) regulates output and error handling. 
0: no output, execution not stopped on error. 
1: only warnings printed, execution not stopped on error. 
2: as 1, but execution stopped on error. 
3: as 2, and output at first initialization, no Minuit output. 
4: as 3, but output at all initialization calls. 
5: full output, i.e. as 4 and MlNUIT output. 

(D=l) regulates information in the event record. To be given as Ii ep ton + 

10 x Ishower, where l\ ev ton = 0/1 inactives/actives the scattered lepton 

(K(i,l)=21/1) and Ishower = 0/1 excludes/includes intermediate partons 

in the parton showers. 

(D=3) choice of frame for the event. 
1: hadronic CM frame, z-axis along exchanged boson. 
2: lepton-nucleon CM frame, z-axis along lepton. 
3: lab system as defined by last user call to LINIT. 
4: as 3, but z-axis along exchanged boson. 

(D=l) regulates the azimuthal angle, 0, of the lepton scattering plane. 
0: no 0-rotation, scattering plane is x — z plane. 

1: random 0-rotation, performed in lepton-nucleon cms for LST(5)> 2. 
(D=l) regulates completeness of Monte Carlo simulation (to speed up pro- 
gram when only partial information is needed). 

1: only kinematic variables generated. 

0: kinematic variables and parton level event generated, optionally including 
QCD effects (cf. LST(8)). Hadronization can be made later by calling 
LUEXEC. 

1: full event generated, i.e. as plus hadronization and decays. 

(D=12) simulation of QCD effects in hadronic final state. 
0: QCD switched off. 

1: first order QCD matrix elements (ME) for gluon radiation and boson-gluon 
fusion. 

2: QCD parton cascade evolution from initial and final quark. 
3: QCD parton cascade evolution from initial quark only. 
4: QCD parton cascade evolution from final quark only. 
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=5: 
=9: 

=12-15: 
Remarks: 



LST(9) : 



= 1 

=2 
=3 
=4 
=5 

=6: 



=9: 



LST(IO) 



=0: 
= 1: 
=2: 
=3: 
=4: 
Remarks: 



LST(ll) : 



LST(12) 
LST(13) 



LST(14) 



=0: 



QCD switched off, but target remnant treatment as in cascade case. 

set by Ariadne [37 when simulating parton emission in the colour dipole 

model fl51] . 

as 2-5, but PS added on event from first order matrix elements. 

See PARL(8),PARL(9) for cut-off on ME. Without PS one may decrease 

the cutoff to maximize QCD emission. 

(D=5) scale in parton showers, i.e. maximum virtuality (mass-squared) of 
parton initiating the shower. Only used when QCD ME are not included, 
i.e. LST(8)=2,3,4. See further section 2.6. 

Q 2 

W 2 

W*Q 

Q 2 (l — x), which is ~ (pj_) from ME at large x 

Q 2 (l — x)max(l,ln-), which is a combination of (p^)-dependence from 
ME at large and small x 

xoW 2 , to represent the mass-squared of the hadronic system whithout the 
(non-perturbative) spectator. A valence-like parton distribution /(xq, Qq) ^ 
(1 — xo) a (a = 4) is used to chose xq representing an original parton. 
PU 4//3 , i.e. similar as in the colour dipole cascade model [37j . 
(D=l) numerical integration of cross section in LINIT. 
not performed. 

50|| , i.e. adaptive Gaussian method. 
, i.e. adaptive Monte Carlo method. 



performed using Gadap 



performed using Riwiad |48" 



performed using Divonne ]4"9" 



performed using Divonne [49 



automatic invocation, 
detailed invocation. 
Integration variables are defined by LST(l), integration region by the 
cuts in array CUT and required accuracy by PARL(15). Result stored 
in PARL(23). Under normal conditions the simpler Gadap routine is ac- 
curate enough. 

(D=0) inclusion of longitudinal structure function, target mass and higher 
twist corrections, see section 2.3, for 7 and j/Z exchange. To be set as 
LQCD+10*LTM+100*LHT where LQCD, LTM and LHT are or 1 to ex- 
clude or include contributions from QCD, Target Mass and Higher Twist, 
respectively. The QCD and target mass parts involve integrals which are 
evaluated at initialization and stored on an x, Q 2 grid used for interpo- 
lation when simulating events. In this mode, the cross section estimates 
PARL(23) and PARL(24) include F L . By setting LQCD or LTM equal 2, 
the corresponding contribution are evaluated by integration for each event 
at the proper x, Q 2 point, which gives a more accurate result. In this 
case, the result is included in PARL(24) but not in PARL(23) since nested 
integrations is too time consuming. 

(D=4) maximum flavour used in sea structure function parametrizations. 
(D=5) heaviest quark flavour allowed in boson-gluon fusion. A threshold 
factor is applied to compensate for neglected quark masses in the matrix 
elements. 

(D=4) treatment of target remnant after removing interacting parton, see 
section 2.7. 

remnant approximated by anti-parton of removed parton, i.e. by q, q, g for 
removed q,q,g. No baryon is produced. 



Lepto 6.5 manual 



=1: 



r(is) 



--2: 

=3: 
=4: 

=0: 



=9 
= 10 
=-4 



for removed valence quark the remnant is a diquark hadronizing into a 
baryon with the Lund model. For a removed gluon, sea quark (g s ), sea 
antiquark the remnant is, respectively, a q v q v q v , qvQvqvQs, <2Wd<2Ws which 
is split into q v q v + q v , q v q v + M(q v q s ), q v + B(q v q v q s ). The (lighter) part 
of the remnant containing one random flavour valence quark q v , takes the 
energy-momentum fraction z given by P(z) = 2(1 — z), i.e. (z) = 1/3. 
as 1, but with P(z) = (a + 1)(1 — z) a with a chosen such that (z) = 
l/(a + 2) = m/(m + M) where m (M) is the mass of the light (heavy) 
remnant subsystem. 

as 2, but using the 'Peterson' function P(z) = N/ (z(l — l/z — c/(l — z)) 2 ) 
with c = (m/M) 2 . 

using LUZDIS, i.e. the fragmentation function chosen in Jetset. 
(D=9) choice of parton distribution functions, xq(x, Q 2 ) and xg(x, Q 2 ), for 
the target hadron, see section 2.4 and LST(16). 

parton density choice and parameters are controlled directly through pa- 
rameters MSTP(51), MSTP(52), MSTP(57), MSTP(58) and PARP(51) in 
common PYPARS in Pythia 5.7 
Eichten-Hinchliffe-Lane-Quigg set 1 



52 



Eichten-Hinchliffe-Lane-Quigg set 2 |52 
Duke-Owens set 1 |53| . 
Duke-Owens set 2 |3|. 
CTEQ2M (best MS fit) H. 
CTEQ2MS (singular at small-x 
CTEQ2MF (flat at small-x 
CTEQ2ML (large A)) fl5~ 
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r(ie) 



CTEQ2L (best leading order fit) 
CTEQ2D (best DIS fit) ||. 

Intrinsic charm quarks only, see section 2.4; to be used with PARL(3)~ 
m c ~ 1 GeV and overall normalisation PARL(12). 
=-5: Intrinsic bottom quarks only, approximated by taking the same distribu- 
tion as for intrinsic charm, but renormalized by m 2 /m 2 ~ 0.1. 
Remarks: Except for intrinsic charm and beauty, the parton densities are obtained 
from subroutine PYSTFU in Pythia 5.7 [p^j ; with parameters set to 
MSTP(51)=LST(15), MSTP(52)=LST(16) and MSTP(58)=LST(12), un- 
less LST(16)=0 when the user directly controls the PYSTFU parameters. 
To evaluate structure functions separately, see subroutines LYSTFU above. 
(D=l) choice of proton parton-distribution-function library, i.e. MSTP(52) 



r(17) : 



r(is) 



in Pythia 5.7 |24|. 

internal Pythia according to LST(15). 

PDFLIB (version 4 or later) with with NGROUP and NSET to be 
given as MSTP(51)=1000*NGROUP+NSET. 

(D=0) regulates varying energies of initial particles from event to event, 
fixed energies as specified in LINIT. 

energies allowed to vary; momenta should be given in P(i,j) with i=l,2 for 

lepton, nucleon and j=1...5 for p x ,P y ,p z , E, m. See also LWBB. Note: this 

option is not fully tested, use ME only with LST(19)=-1. 

(D=2) running of electromagnetic coupling a and choice of W and Z 

masses. 
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LST(19) 



LST(20) 



LST(21) 
LST(22) 
LST(23) 



LST(24) 

LST(25) 
LST(26) 
LST(27) 



=4: 
--5: 

=6: 



=0: a fixed to value at Q 2 = given in PARL(16), Z and W masses given 

independently by PMAS(23), PMAS(24) in /LUDAT2/. 
= 1: as 0, but W, Z masses calculated from standard model using sin 9w, ctjGp 
and radiative corrections; i.e. PARL(5), PARL(16), PARL(17), PARL(18) 
=2: as 1, but a varying with Q 2 using ULALEM in see MSTU(lOl). 
: (D=-10) regulates use of grid to store probabilities for qg- and gg-events. 
=-10: the probabilites are obtained from an automatic grid, but close to the 

saturation limit they are recalculated to avoid interpolation errors. 
=-1: no grid, necessary integrals calculated for each event, more time consuming 

but higher precision. 
=0: user defined grid to be read in free format, see comments in subroutine 
LWEITS. 

= 1: grid suitable for lepton beam energy < 300 GeV on fixed target. 
=2: grid suitable for lepton beam energy < 1000 GeV on fixed target. 
=3: grid suitable for ep collisions in HERA. 
=4: grid suitable for ep collisions in LEP+LHC 

=10: grid chosen automatically after the kinematic region specified by user. 
: (D=5) scheme for cut-offs against divergences in the QCD matrix elements 

f3~4| , cf. section 2.5 and PARL(8),PARL(9). (to^ is the invariant mass of 

any parton pair). 

W 2 (or JADE) scheme, to|- > y cu tW 2 
Q 2 scheme, to?- > y cut Q 2 

mixed scheme, i.e. m 2 - > C\Q 2 for partons i,j that are not spectators and 
2pi ■ P > C2Q 2 /x = C2(W 2 + Q 2 ) else. Virtuality scale for matched parton 
showers in g-event is C\Q 2 . 

as 3, but with virtuality C2Q 2 jx for initial state parton showers. 



■s scheme, 



s > s r 



and Zq <^ Zq y rn%i% 



< 1 



To regulate the divergence, 



is changed, see PARL(27). 



as 5, but Smin is changed to regulate the divergence, see PARL(27). 
error flag, =0 for properly generated event. Nonzero value indicates in- 
correct event that should be rejected; may occur in case of user supplied 
kinematic variables or variable beam energies, cf. LST(2) and LST(17). 
specifies chosen target nucleon in current event. 
=1: proton. 
-2: neutron. 

specifies process simulated. 
=1: electromagnetic (EM), i.e. 7 exchange. 
-2: weak charged current (CC), i.e. W ± exchange. 
=3: weak neutral current, i.e. Z° exchange. 
=4: neutral current (NC), i.e. 7/Z exchange. 

specifies first order QCD process in current event. 
= 1: g-event, i.e. no first order QCD. 
-2: gg-event, i.e. gluon radiation in first order QCD. 
=3: gg-event, i.e. boson-gluon fusion in first order QCD. 

specifies flavour of struck quark in current event: 1—d, 2=u, 3=s, 4=c, 
5=6, -l=d, -2=u, -3=s, -4=c, -5=6. 

entry line in event record of outgoing struck quark. In parton shower case, 

quark at boson vertex before final state shower. 

signals split of non-trivial nucleon remnant, cf. LST(14). 
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LST(28) 

LST(29) 
LST(30) 



=0: no split, simple diquark or LST(14)=0. 

= 1: split into parton and particle, qq + M or q + B, occurs when sea (anti)quark 

removed through the interaction. 
=2: split into quark and diquark, q + qq, occurs when a gluon is removed. 

specifies the frame in which the current event is given with code as for 

LST(5), cf. LFRAME. 

specifies azimuthal angle rotation with code as for LST(6), cf. LFRAME. 
specifies chosen helicity of beam lepton in current event. 
-1: left-handed. 



+1: 



right-handed. 

specifies chosen variables for the simulation/integration. 
x and Q 2 . 
x and y. 
x and W 2 . 

Internal flag, or 1 for simulation and integration, respectively. 
Reserved for internal test. 

(D=l) switch for soft colour interactions, 0=off, l=on; see section 2.8. 
(D=l) switch for new sea quark treatment, 0=off, l=on, see section 2.7. 
(Inactive for scattering on intrinsic charm.) 
LST(36)-LST(40) : unused at present. 



LST(31) 



LST(32) 
LST(33) 
LST(34) 
LST(35) 



PARL(l) 
PARL(2) 
PARL(3) 

PARL(4) 

PARL(5) 
PARL(6) 



PARL(7) 



(D=l.) number of nucleons in target nucleus, i.e. A. 
(D=l.) number of protons in target nucleus, i.e. Z. 

(D=0.44 GeV) width of Gaussian distribution for the primordial transverse 
momentum k± of partons in the nucleon. 

(D=0.75) probability that a wrf-diquark in the target remnant has spin and 
isospin equal zero, i.e. I=S=0. 
(D=0.2319) sm 2 6 w (Weinberg angle) g|. 

(D=0.) polarization of lepton beam; should be set to 1 — 2Pl = 2Pr — 1 
where Pl,r is the probability of left(right)-handed helicity, i.e. — 1(+1) 
corresponds to a fully left(right)-handed polarized beam. 
(D=0.5) probability for soft colour interactions between parton pairs, see 
section 2.8. 

PARL(8),PARL(9) : (D=0.04,4) cut-offs against divergences in the QCD matrix elements, 
cf. section 2.5 and LST(20). Suitable values are given in parenthesis. 



LST(20) = 
LST(20) = 
LST(20> 



=1,2 
=3,4 
=5,6 



PARL(8> 
PARL(8> 
PARL(8> 



y cut and PARL(9)=minm ii in GeV. (0.005,2) 
c 2 and PARL(9)= c x . 



and PARL(9> 



(0.05,2) 
s min in GeV 2 . (0.04,4) 



Remarks: These are starting values when integrating first order QCD matrix ele- 
ments, but the effective cut used, PARL(27), is automatically increased in 
order that the QCD probabilities do not exceed unity. 
Not used at present. 

(D=0.01) required relative accuracy for one-dimensional integration, used 
for first order QCD matrix element weights and longitudinal structure 
function integrals. 

(D=0.01) probability f3 2 for an intrinsic charm quark-antiquark pair in the 
proton (cf. section 2.4). 

(D=0.1) internal parameters used for adjustment of y cu t for integration of 
QCD matrix elements. 



PARL(10) : 
PARL(ll) : 



PARL(12) : 
PARL(13) : 
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PARL(14) 

PARL(15) 

PARL(16) 
PARL(17) 
PARL(18) 
PARL(19) 
PARL(20) 

PARL(21) : 

PARL(22) : 
PARL(23) : 



(D=0.35 GeV) width of Gaussian distribution in transverse momentum 
when a non-trivial target remnant is split into two particles, cf. LST(27). 
(D=0.01) required relative accuracy for two-dimensional integration to get 
total cross section, cf. LST(IO). 

(D= 7.29735 x KT 3 ) finestructure constant a [0, cf. LST(18). 
(D= 1.16639 x KT 5 GeV -2 ) weak Fermi constant G F 0. 
(D=0.044) Ar from radiative corrections [I9|. 

(D=0.03) scale k 2 in GeV 2 for higher twist correction | 22| , see eq. (f2~7j). 



(D=0.1 GeV) in case of a more complicated nucleon remnant than a simple 
diquark, the constituent masses of quarks and diquarks in the remnant are 
decreased with PARL(20) and 2*PARL(20) respectively. 
IP ■ k, where P is proton and k lepton four-vectors, equals invariant mass 
squared, s, when masses are neglected. 
IP ■ q, where P is proton and q is boson four-vectors, 
cross section in pb, corresponding to the kinematic region allowed by the 
cuts in array CUT, obtained by numerical integration in LINIT call, cf. 
LST(IO). If LQCD=2 or LTM=2, cf. LST(ll), the corresponding contri- 
bution to F L is not included. 
PARL(24) : Monte Carlo estimate of the cross section in pb associated with the gen- 
erated event sample, taking CUT into account. Set to zero by LINIT call 
and updated with each event generated, hence accuracy improves with 
statistics as 1/ vN and only final value should be used. 
PARL(25) : value of a s in current event. 

PARL(26) : value of A obtained from structure functions, which is used in a s and in 
the initial state parton shower. 

PARL(27) : depending on LST(20), present value of the y cu t~, z q - or s-cut for first order 
QCD. These are given by PARL(8) and PARL(9), but modified internally 
to prevent QCD weights to exceed unity. 

PARL(28),PARL(29),PARL(30) : values of x p , z q and cf) in first order massless QCD ma- 
trix elements, section 2.5, for current event if it is a qg- or gg-event, see 
LST(24). For a q— event they are set to 1.0, 1.0 and 0.0. 



X : Bjorken-x, i.e. x = Q 2 /2P ■ q. 

Y : standard y variable, i.e. y = P ■ qj P ■ p e . 

W2 : mass-square of hadronic system, i.e. W 2 = (P + q) 2 . 

Q2 : momentum transfer squared, i.e. Q 2 = —q 2 = —{p e — pe) 2 - 

U : energy transfer variable v = P ■ q/m p . 

Remark: for details see section 2.1. 

COMMON /LFLMIX/ CABIBO(4,4) 

Purpose: Contains the Cabbibo-Kobayashi-Maskawa matrix elements squared for 
flavour mixing. First index corresponds to the up-quark flavours u, c, t, h 
and second index to the down-quark flavours d, s, b, I. The default values 
are CABIBO/.95,.05,2*0.,. 05, .948,.002,2*0., .002,. 998,4*0., 1./ 

COMMON /LOPTIM/ OPTX(4),OPTY(4),OPTQ2(4),OPTW2(4),COMFAC 

Purpose: parameters to optimize simulation efficiency, see section 2.2. Default values 
are good for normal usage, but improvements may be possible under, e.g., 
unusual kinematic conditions. For changes, see their detailed meaning 
given by comments in the code of subroutine LEPTOX. 
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COMMON /LBOOST/ DBETA(2,3),STHETA(2),SPHI(2),PB(5),PHIR 

Purpose: rotation and boost parameters for transforming event between frames. 

DBETA(ij), STHETA(i) and SPHI(i) are the boosts, polar and azimuthal 
angles for transforming from lepton-nucleon cms to lab (i=l) and from 
hadronic cms to lepton-nucleon cms (i=2). For these cases, the order of 
the transformations should be first rotation in 9, then in followed by 
the boost where j— 1,2,3 corresponds to x,y,z components. Should only 
be used via subroutine LFRAME to set flags properly, cf. LST(5) and 
LST(6). 

COMMON /LGRID/ NXX,NWW,XX(31),WW(21),PQG(31,21,3),PQQB(31,21,2), 
& QGMAX(31,21,3),QQBMAX(31,21,2),YCUT(31,21),XTOT(31,21),NP 

Purpose: information for simulating first order QCD processes. 

Remarks: Probabilities PQG and PQQB for gluon radiation and boson-gluon fusion, 
i.e. qg- and gg-events, are stored on a grid in x, W (or x, y for LST(19) = 10) 
with NXX, NWW points defined by XX and WW. Indices are for x, W (y) 
and helicity contributions. The cut on parton pair masses y^ is also stored 
and maxima for the Monte Carlo rejection technique. The grid content 
is set up by LWEITS and can be printed by LPRWTS. See LST(19) and 
section 2.5 for physics and methods. 

COMMON /FLINFO/ RFLQ,RFLG,RFLM,RFLT 

Purpose: gives the relative contributions from the different parts (quark, gluon, mass, 
higher twist) of the longitudinal structure funtion to the differential cross 
section at the x, Q 2 point of the current event, cf. section 2.3. 

COMMON /LYPARA/ IPY(80),PYPAR(80),PYVAR(80) 

Purpose: parameters for parton cascade routines. This is common PYPARA from 



Pythia 4.8 and described in detail in |[40|| . Only those parameters used 
with Lepto in parton cascade mode are commented here (default values 
as in Pythia 4.8 if not given). 

Parameters: 

IPY(8) : set to LST(12) in UNIT. 

IPY(13) : set to zero in LINIT if LST(8)=3,5, 13,15. 

IPY(14) : set to zero in LINIT if LST(8)=4,5, 14,15. 

IPY(11),IPY(15),IPY(34),IPY(40)-IPY(42),IPY(47),IPY(48) : are used. 

PYPAR(8),PYPAR(11)-PYPAR(16) : are used. 

PYPAR(21) : (D=PARL(26)) A QCD in initial state parton shower. 

PYPAR(22) : (D=l. GeV 2 ) cutoff for initial state parton shower. 

PYPAR(23),PYPAR(24) : are used. 

PYPAR(25),PYPAR(26) : (D=2*l.) multiplies the chosen scale, cf. LST(9), to give max- 
imum virtuality for final and initial state showers, respectively. 

PYPAR(27) : (D=l.) multiplies virtuality Q 2 for a s and structure functions in initial 
state showers. 

PYVAR(1)-PYVAR(5) : are used. 

COMMON /LUJETS/ N,K(4000,5),P(4000,5),V(4000,5) 

COMMON /LUDAT1/ MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
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Purpose: 



LUJETS contains the record of the generated event and is essential for us- 
ing the results. LUDAT1 contains switches and parameters that are, e.g., 
essential to control final state parton showers, a s evaluation and hadroniza- 



tion. These common blocks are described in the Jetset manual [24]. 



Short description of all common blocks: 



Purpose 

main common block for user control of program, see above, 
quark flavour mixing parameters (KM- matrix), see above, 
parameters to optimize simulation efficiency, see above, 
boost and rotation parameters, see above, 
grid for first order QCD event simulation, see above, 
internal parameters and variables; charges, couplings, cross section weights, 
internally used; basic system in different frames, effective limits on kine- 
matic variables, 
internal for output control, 
internal counters for integrand evaluations, 
grid with integrals for longitudinal structure function, 
relative size of longitudinal structure function contributions, see above. 



Common 
LEPTOU 
LFLMIX 
LOPTIM 
LBOOST 
LGRID 
LINTER 
LINTRL 

LPFLAG 
LINTEG 
FLGRID 
FLINFO 

LYPARA parameters for parton cascade evolution [40], see above. 
LYPROC,LYINTl internal for parton showers [4T1 ]. 

LMINUI,LMINUC input parameters and character names for Minuit |20 
LM.... internally used in adaptation of Minuit routines. 

GAD API internal in Gadap integration routines. 
PARAMS , ANSWER input /output for Riwiad integration @. 
STORE, STORE1, OPTION, RANDOM, INTERN internal for Riwiad |S 



BNDLMT, SAMPLE, PRINT input for Divonne integration |49 
ARDAT1 



parameters in the Ariadne MC for dipole radiation [47 



LUDAT1,LUDAT2 switches, parameters, particle data in Jetset |24 . 
PYPARS switches, parameters used for parton densities in Pythia 5.7 [2~3|. 
LUJETS contains generated event, see |24 
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4 Usage and availability 

Lepto 6.5 should be loaded together with Jetset 7.4 and Pythia 5.7 p4j (the latter 
to access parton density parametrizations) . The ordinary gamma function, GAMMA(X), 
is called and must be supplied (usually available in FORTRAN77). Access to the CERN 
library is not necessary, but gives access to the Pdflib, Riwiad and Divonne program 
packages as well as the TIMEX routine. The program is a slave system, which the user 
must call from his own steering program. 

Information about the program, its update history, source code, example jobs etc. can 
be obtained on request from the authors or directly via the WWW on the Lepto home 
page, tittp : //www3 . tsl ,uu. se/thep/lepto/ . 
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